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Abstract 

This work investigates the first correction to the equilibrium phase space distribution and its 
effects on spectra and elliptic flow in heavy ion collisions. We show that the departure from 
equilibrium on the freezeout surface is the largest part of the viscous corrections to V2{pt)- However, 
the momentum dependence of the departure from equilibrium is not known a priori, and it is 
probably not proportional to p'^ as has been assumed in hydrodynamic simulations. At high 
momentum in weakly coupled plasmas it is determined by the rate of radiative energy loss and 
is proportional to pj, . The weaker px dependence leads to straighter V2{pt) curves at the same 
value of viscosity. Further, the departure from equilibrium is generally species dependent. A 
species dependent equilibration rate, with baryons equilibrating faster than mesons, can explain 
"constituent quark scaling" without invoking coalescence models. 
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I. INTRODUCTION 



When two ultra-relativistic nuclei collide, they leave behind a region of high energy- 
density QCD matter, whose properties we would like to understand better. The initial 
geometry of the QCD matter is set by the overlap region of the two colliding nuclei. Gen- 
erally, the nuclei collide at finite impact parameter rather than head-on. In this case the 
initial geometry is not a disk, but is an "almond shaped" ellipse. (The short and long axis 
of the initial almond are taken as the x and y axes respectively.) The production mechanism 
of the QCD matter is local and knows nothing of this global geometry. Therefore, to a first 
approximation the initial stress tensor will be locally azimuthally symmetric. Subsequently, 
if there are no reinteractions the produced matter will free stream to the detector; the ini- 
tial geometry will have no influence on the evolution, and the angular distribution of the 
final observed hadrons will also be azimuthally symmetric. On the other hand, if there are 
strong interactions which maintain local thermal equilibrium, the pressure gradients in the 
X direction will be larger than in the y direction, an anisotropy in the collective flow will 
develop, and ultimately an anisotropy in the momentum spectrum of the final hadrons will 
be observed. 

The final momentum anisotropy is characterized experimentally by f2, the second har- 
monic of the azimuthal distribution of the produced particles with respect to the reaction 
plane. Experimentalists have measured V2 as a function of transverse momentum px-, par- 
ticle type, and impact parameter These results are surprisingly well described by 
ideal hydrodynamics |5i], which amounts to the approximation that the interactions are fast 
enough to maintain the matter in equilibrium from an early time until hadronic freeze-out. 
There are some limits to this success. First, the measured V2 falls below the ideal hydro- 
dynamic prediction for momenta larger than pt ~ 2.0 GeV. Second, the hydro fit fails to 
reproduce certain relative trends observed in the baryon and meson elliptic flows. These 
trends are compactly summarized by "constituent quark scaling" which generally has 
been attributed to a kind of coalescence of constituent quarks UEJ. Here we will argue 
that the first corrections to equilibrium can clarify both of these shortcomings without the 
need for a coalescence model. 

To quantify the corrections to ideal hydrodynamics it is important to study nonideal 
viscous) hydrodynamics. In the last two years there has been a major push in this direction 
13H20|. These studies have used various formalisms and have studied variations of V2{pt) 
with respect to the input shear viscosity, the model for the initial geometry, and various 
other nuisance parameters. However, we want to point out here that these studies have 
all made a common assumption about the way that the asymmetry in the stress tensor is 
manifested in the particle distribution after freezeout. In particular, the particle distribution 
after freezeout is locally of the form / = /o + 5/ where /o is the equilibrium distribution and 
5f is the first correction. All groups have assumed that 5f{p) oc p^/o and that the coefficient 
of proportionality is independent of particle type. 

In this paper we will argue that this assumption matters, and that it is far from secure. 
After an overview of the issue in the next section, in Section UlI] we will discuss the physics 
which establishes the momentum dependence of Sf and its behavior in several theories. We 
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will see that while the most studied theories give Sf oc p^fo, the most QCD-like theories do 
not. Then we explore the behavior of multi-component plasmas in Section HVl We see there 
that the viscous corrections Sf{p) for different species are generically different. This fact can 
account for the "constituent quark scaling" observed in the baryon and meson elliptic flows 
without any reference to the hadronization process. We then make our concluding remarks. 
Some technical material is postponed to appendices. 

Throughout, we will denote 4- vectors with capital letters P, Q and use p, q for their 3- 
vector components, Ep,Eq for their energy components, and p,q for |q|. Our metric 
convention is [-,+,+,+], so that m^m^ = —1. We use tilde to indicate momenta scaled by 
temperature, p = p/T. We will mostly write Hp for the equilibrium distribution function. 
Up = l/(exp(p/T) =F 1) but will occasionally use fo{p) when common convention dictates its 
use. The appropriate statistics will be clear from context. 



II. OVERVIEW 



The energy momentum tensor is given by the sum of its ideal and dissipative parts^ 

T^"" = (e + V)u''u'' + Vg^"" + 7^'^^ (1) 
and obeys the equation of motion. 



df,T>"' = 



(2) 



In the first-order (or Navier-Stokes) approximation the dissipative part of the stress energy 
tensor in the local rest frame is 



(3) 



where rj is the shear viscosity, and we use (. . .) to indicate that the bracketed tensor should be 
symmetrized and made traceless. It is well known that the first order theory is plagued with 



difficulties such as causality violations and instabilities 2m |22| . In order to circumvent these 



issues a second order theory is required. The most commonly used second order relativistic 



viscous hydrodynamics is due to Israel and Stewart |23|. For technical reasons we use 



a theory developed by Ottinger and Grmela 2J, |25|. The two theories are qualitatively 
the same {i.e. for sufficiently small relaxation times they both approach the first order 
theory). To streamline the presentation we postpone the details of our hydrodynamic model 
to Appendix lA] and refer to previous work [l7|. 

The solutions to the hydrodynamic equations yield the underlying temperature and flow 
profiles in the presence of viscosity. Particle spectra are then computed using the Cooper- 
Frye f26| formula 



E 



V 



d^p (27r) 



f{p)p^da^ 



(4) 



^ We use Landau-Lifshitz conventions to fix e,u^ in terms of four components of T^". The other six 
independent components of T'"^ can always be accommodated by a tt'"^ satisfying u^tt^"^ = 0. 



3 



where p = p/T and is the freeze-out hypersurface taken as a surface of constant energy 
density in this work. For a system out of equihbrium f{p) is not the equihbrium distribution 
function but also contains viscous corrections, 

f{p) = fo{p) + Sf{p), (5) 

where /o is the ideal Bose/Fermi distribution function. The form of 6f is constrained by 
the requirement that T*-' be continuous across the freeze-out hypersurface: 

Dropping 5f from the final particle spectra is inconsistent as it leads to a discontinuity in 
T^'^ . The form for 5f which satisfies continuity in the local rest frame is proportional to 
P^P'T^ij and is traditionally parametrized by xip) ^ 

5/(p) = -n,{l±n,)x{p)^ (7) 
= -X{.p)np{l ± np)py {diUj) , (8) 

where we have distinguished x{p) and xip) = x{p)p^P' {diUj) by the argument of the func- 
tion. One moment of x{p) is fixed by the shear viscosity (see below) but otherwise x{p) 
is an arbitrary function of p. To date all works on viscous hydrodynamics have taken the 
quadratic Ansatz and have usually worked in a Boltzmann approximation 

xip) oc . (9) 

Let us look at typical results for V2{pt) as shown in Fig.[Tl The curve labeled 'Ideal' shows 
the result using ideal hydrodynamics {i.e. rj/s = 10'^). The curve labeled '/o' shows the 
resulting elhptic flow from the viscous evolution (the solution of Eqs. (I2]),(l3]) and (IA2I) ) but 
without including the viscous correction to the distribution function. In other words, this 
shows how the viscous correction to the temperature and flow profiles manifests itself in the 
particle spectra. Only modest corrections to the spectra are found. As already emphasized, 
this result is unphysical since dropping 6f violates continuity of the stress tensor. Last, 
the curve labeled '/o + 5/' also takes 6f into account, using the quadratic Ansatz. The 
viscous correction to the distribution function dominates the reduction in V2 at large pt. 
That means that the 6f term is responsible for a significant part of the effects of viscosity 
in the particle spectra. 

This being the case, it is imperative to perform a systematic study on the form of the 
viscous correction as well as its effect on elliptic flow. Most of this paper will discuss the 
form of the viscous correction appearing in weakly-coupled QCD. Although this is not a 
theory of hadronizing QCD, it is one theory where quantitative first-principle calculations 
can be performed. One of our major findings is that not all models of energy loss give the 
same predictions for the off-equilibrium distribution function. 



^ In actual simulations tt*-' is treated as a dynamical variable in a second order fluid formalism. Then to first 
order one can make the replacement, (diUj) — >■ —-Kij/irj. There has been no attempt to systematically 
include 5f through second order in hydrodynamic simulations. 
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Pt [GeV] 



FIG. 1: Typical results for V2{pt) from a viscous hydrodynamic model employing the quadratic 
Ansatz. The run parameters are rj/s = 0.08, Tfrzom = 140 MeV and p = l/3e. Further details are 
in Appendix lAl 



III. FORM OF 5f IN SEVERAL THEORIES 



In this section we consider a number of theories, to show that while the dependence 
xip) c<: is expected in some cases, other functional dependence is expected in others, 
including weakly coupled QCD and a hadron (resonance) gas. The theories where we can 
make a definite statement about the functional form of 5f are all described by kinetic theory. 
Since freeze-out is defined as the point where scatterings go from being common to being 
rare on the time scale of the evolution of the system, we generally expect that, just before 
freezeout, kinetic theory should be a reasonable description. 

Within kinetic theory, the distribution function /(p, x) is determined by a Boltzmann 
equation, 

{dt + v^-d^)f{p,x) = -G[f,p], (10) 
where C[/, p] is the collision operator. In equilibrium the distribution function obeys 

<p^^) = ^.p^u^^tXnt,-)^! ^ cKp]=o. (11) 

To determine the first viscous correction 5/, we work in a vicinity of the local rest frame 
= (1, u^{x, t)), and substitute / = n{p, x) + 6f into Eq. (|TOl) keeping terms first order in 
the spatial derivatives 

^n,(l ± n,) {d,u,) = -C[Sf,p] . (12) 

Here C[6f,p] denotes the linearized collision operator, i.e. the collision operator expanded to 
first order in 6f. In writing Eq. f|T2|) we have used ideal hydrodynamics and thermodynamic 
relations to rewrite time derivatives as spatial derivatives, and we have neglected gradients 
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proportional diU^ which are responsible for the bulk viscosity 27|. Eq. ( !T2|) is an integral 
equation for 5f which can be solved by various methods. 

Since the first viscous correction is a scalar and must be proportional to the the strains, 
the most general form for the viscous correction in the local rest frame can be parametrized 
by the function xip) as in Eq. ([71). Close to equilibrium the first viscous correction 5f 
determines the strains 



TT 



-2.,(9V> = y p^i^^/. (13) 

which ultimately yields a relation between between the shear viscosity and the viscous 
correction x{p) 

This is the only general constraint on the functional form of the viscous distribution function. 
To proceed further we must specify completely the form of the linearized collision operator, 
which we will do in the context of various model theories. 



A. Simplest model: relaxation time approximation 

The simplest model (really a cartoon) for the collision operator is the relaxation time 
approximation, 

^df,p\ = — , . , (15) 

where th is the momentum dependent relaxation time to be specified. Substituting this 
form for the collision operator into Eq. f lTS]) . and working in a Boltzmann approximation 
np{l ± Up) — )■ Hp yields the following form for 6f: 

(5/ = -^^n,pV(9,«,). (16) 

Note however that the relaxation time is in general energy dependent. In different theo- 
ries, Tji{Ep) might show different functional dependence on Ep. Without details about the 
dynamics of the theory in question, we can only parametrize the viscous correction. Here 
we will discuss a massless classical gas where Up = e~^/^ and parameterize the relaxation 
time (or the distribution function) with a simple power law 

Sfip) = -npXiPWP^ {diUj) , 

xiP) = C{a)p'-^ . (17) 



The constant, C(a), is determined through Eq. (IT4] 
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FIG. 2: Left: V2{pt) using the Linear or Quadratic Ansdtze for the distribution function. Right: In- 
tegrated V2 versus centrahty showing independence from the precise form of the viscous correction. 
Run parameters can be found in Fig. [H 

There are two limiting cases for the functional form of the the relaxation time approx- 
imation, a = and a = 1. The momentum dependence of the relaxation time in these 
extreme cases is 

, , J p a = (quadratic ansatz), , . 

tr p oc <y ^^^g^ a = 1 (linear ansatz). 

Most theories will lie between these two extreme limits.^ Loosely speaking, if the energy loss 
of high momentum particles grows linearly with momentum, ^ oc p one expects a relaxation 
time independent of momentum, tr oc p^. On the other hand if the energy loss approaches 
a constant ^ oc const, the relaxation time will grow with the particle momentum tr oc p. 

Fig. [2] shows the elliptic flow computed using these two functional forms for the first 
viscous correction. It is important to emphasize that shear viscosity is the same in both 
cases. Examining these figures, we see that the integrated elliptic flow is largely insensitive 
to the functional form of the first viscous correction. This is because the integrated V2 is 
primarily determined by the hydrodynamic variables e, u^^'k^'^ which are independent of the 



functional dependence of the relaxation time 27|. The differential elliptic flow V2{pt) is 



sensitive to the rate of equilibration especially above pr — 1-2 GeV. 



B. Scalar \<j)'^ theory 



Scalar field theory has been described at length by Jeon 28|], who rigorously derived the 



Boltzmann equation and its collision kernel and then solved for x{p) numerically. But if we 



^ There are exceptions to this rule. For instance, in a gas of Goldstone bosons far below the symmetry 
breaking scale one expects a — 2, since the cross section grows rapidly with energy, cr ~ S^/A"'. 
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make the approximation of Boltzmann statistics, we can actually solve for Sf in closed form. 
First the non-linear Boltzmann equation with Bose-Einstein statistics is 



C[f,p] = / Tpk^p,k' [fpfkil + + fk') - fp'fk'il + + fk)] 

Jk,p',k' 

where the transition rate (including a final state symmetry factor) is 



(20) 



1 IMI' 



and we have used the traditional short hand, jp = J (f^- After linearizing, f{p) — )■ 
Up + np{l ± np)x{p) the linearized collision integral is 



p'-pi 

(diUj) 



TEp 



/ Tpfc^p/fc/ UpUkil + np,){l + Uk') [x{p) + x{k) - x{p) - x{k')] . (22) 

Jk,p',k' 



At this point we will make the Boltzmann assumption by neglecting the stimulation factors, 
(1 + Up) — 7- 1, and using, Up = e~^^'^. Then we will try a solution of the form 

X{p) = Cpy {d.u,) /T' , (23) 

i.e. assuming that x(p) oc or a = 0. Substituting this form into the integral equation 
(Eq. (IT^ and Eq. ( 12^ ). and performing the integrals yields (see Appendix iBj) 

{d.u,) ^ = C (d^u,) • (24) 

Thus, taking C = 3847r^/A^, the quadratic form in Eq. f l2^ has provided an exact solution 
solution to the linearized integral equation. The viscosity is rj = 15367rT^/A^ in a Boltzmann 
approximation. 

Physically, this happens because of the form of the scattering cross- sect ion. Since a oc 
A^/ s and s oc p, the cross-section scales as the inverse of the particle's energy. The typical 
scattering is nearly randomizing, but high energy particles undergo fewer scatterings than 
low-energy ones. Therefore we find the same functional form as for momentum diffusion but 
for very different reasons. 

This example from scalar field theory, together with the example of momentum diffusion 
(see Section [III C II below), are the reason that most people assume the quadratic ansatz, 
x{p) oc should hold. 



C. Weakly coupled pure-glue QCD 

In this section we will use the Boltzmann equation for pure-glue QCD in three approx- 
imation schemes to calculate the first viscous correction. First we will consider a leading 
log(T/m£)) approximation where the dynamics can be summarized by a Fokker-Planck equa- 
tion which describes the momentum diffusion of quasi-particles. In this limit we will find 
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that the viscous correction is quadratic at large momentum, x{p) oc p^. Next we will con- 
sider the QCD Boltzmann equation but consider only 2 — )• 2 collisions and neglect collinear 
radiation. In this limit, we will find that the viscous correction at large momentum behaves 
x{p) oc P^/ log(p). Finally, we will also include collinear radiation in the Boltzmann equa- 



tion as is necessary in a complete leading order treatment 29|, |30|] . We will find that collinear 



radiation controls the relaxation of the high momentum modes and asymptotically we have 
x{p) oc p'^/^, where the coefficient of proportionality is set by the rate of transverse momen- 
tum broadening, q. The impatient reader may skip to Section Till C 41 which summarizes the 
results of these three approximation schemes. 



1. Momentum diffusion in a leading log treatment 

In a leading log approximation, \og{T /m£,) is considered a large number and the dynamics 
describes soft Coulomb scattering. Each soft collision involves a small momentum transfer of 
order q ~ gT, but these collisions happen relatively frequently at a rate of ~ g^T (neglecting 
logarithms). Thus a typical particle with momentum T will diffuse in momentum space and 
equilibrate on a time scale of ~ g'^T . The resulting Boltzmann equation linearized around 



equilibrium can be written as a Fokker-Planck equation [31|, 132 



dt6f + Vp- djf = Tfi 



_d_ 

dp'' 



np{l + Up) 



+ gain terms , 



(25) 



where n is the drag coefficient of a high momentum gluon in this approximation scheme 33 
3 



dp 



with 



/i 



247r 



-^T^log 



T 

mo 



(26) 



The precise form of the gain terms has been given in |3ll . l32| . but only involves the i = 
0,1 spherical harmonic components of Sf{p), i.e. f dQp6f{p) and / dflpp6f{p) . In 
the hydrodynamic limit considered here Sf{p) is proportional to a traceless rank 2 tensor 
(^pipj _ (5«jy3) and these gain terms vanish. Substituting the form of Eq. ([7]) into Eq. ( 125|) 
leads to the following equation for x{p)- 



np{l + np)|; = TjjHpil + Up 



d^ 



dp"" 



+ 



l + 2np 2W ^ 6 I ^^^^ 
T p J dp p^ 



We are not aware of a closed form solution to this equation, but we can find a solution 
for x{p) large momentum. Making the approximation 1 + 2np ^ 1, we find that 



x{p) 



p" 



2T/i 



(28) 



solves this equation. This is the well known quadratic Ansatz. 
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2. Boltzmann equation with 2 — )• 2 collisions 



We next will consider the QCD Boltzmann equation but we will neglect coUinear radi- 
ation. We emphasize that this is not a consistent approximation scheme. Nevertheless, it 
illustrates clearly the relative roles of hard collisions and inelastic processes in determining 
the functional form of x{p) i^i the relevant sub-asymptotic regime. 

The linearized Boltzmann equation is the same as Eq. fl22l) . but the squared matrix 
element is 

|A<|= = 8/Cl(3-|-f-|), (29) 

which describes 2 — )■ 2 gluon scattering after summing over all spins and colors and dividing 
by the gluon degeneracy factor 2d a- These matrix elements must be dynamically screened 
using Hard Thermal Loops. A procedure which is consistent at leading order (where the 
Debye mass is small) but which makes a reasonable estimate when the Debye mass is not 



small has also been described in [35|, and we can follow exactly the numerical procedure of 
that reference to find x(p)-^ We can also study the asymptotic behavior more directly. At 
asymptotically large momentum where log(p) may be considered large, Appendix |B] shows 
that 

X{P) OC :j — —j^ . (30) 

log(p/T) 

The constant in front of the log is related to {dE/dt)^, the rate of collisional energy loss of 
a gluon with momentum p, 

2 

^^^^ = 2twm; • ^^'^ 



In a leading \n{p/T) approximation the loss rate is [34, l36 



^\ =^TMogf^), (32) 
dt /p 487r ^ \tJ ' ^ ^ 

as is rederived in Appendix |Bl The above asymptotic form agrees well with the numerical 
solution of the Boltzmann equation. 



3. A leading order treatment at asymptotically large momenta 



Earl y c alculations of the shear viscosity in pure-glue QCD found x{p) ^ V^i that is. 



a 



[371, l38|. However this is because they were leading-log treatments, which reduced 



29, 30 that inelastic number 



to momentum diffusion discussed above. It was realized in 
changing processes are only suppressed by a log, but are enhanced at large energy E hy & 
factor of (E/TY^'^ and dominate equilibration for E/T > log(l/5'). 

This should not be a surprise. After all, if we think about "equilibration" (energy loss) in 
QED, we find that although the leading order mechanism for the energy loss of a high energy 



Some minor technical difficulties are discussed in the next section. 
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electron is ionization (elastic scattering), bremsstrahlung actually dominates the loss rate. 
This is the case because in bremsstrahlung the energy lost per scattering can scale with the 
incident energy, rather than being incident energy independent as is the case with ionization. 
As a result, the penetration depth of an electromagnetic shower scales only logarithmically 
with the incident energy, i.e. the relaxation time is constant up to logs, oc E^. If the 
same behavior occurred in QCD we would expect the linear Ansatz to hold, a = 1. 

The current understanding of energy loss in perturbative QCD is that the high energy 
behavior lies between these extremes. High-energy particles in a QCD plasma lose energy 
predominantly by inelastic gluon radiation and the time scale for energy loss is short com- 
pared to the time scale for momentum diffusion ("jet broadening"). In particular it was 
shown by Baier et al that for E the rate of (inelastic) energy loss scales with the inci- 
dent energy as dE/dt oc i?^/^, with the half-integer power arising from the LPM suppression 



39I . |40|. This implies a "relaxation time" which scales as r/j ~ E/ [dE/dt) oc i?^/^, and 
therefore a = 1/2 3^. Let us see how this emerges in the behavior of pure-glue QCD. 

The point is that the Boltzmann equation for a gluon plasma possesses both an elastic 
scattering term and an inelastic effective 1 — )■ 2 scattering term, 

dtf + ■ d^f = -C'^V] - C'^'V] ■ (33) 

This equation was first solved at leading order in by Arnold, Moore and Yaffe to determine 
the shear viscosity [sSj]. Their approach involved writing a multi-parameter Ansatz for xiP) 
in terms of a basis of test functions. While the determination of rj improves quadratically 
with the test function basis, the determination of xiP) improves only linearly. Therefore to 
get good accuracy out to p = 15T requires the use of a large basis of functions. We find a 
basis of eight functions is sufficient and the p^/^ behavior is already clear with such a basis. ^ 
We can also directly establish the asymptotic form of the solution. At asymptotically 
high momentum near collinear bremsstrahlung dominates the equilibration of gluons. We 
therefore look at the Boltzmann equation including only 1 — )■ 2 splittings, 

dJ + v^-d.f = -C'^'[f]. (34) 

The relevant collision integral for near collinear joining and splitting of gluons at leading 
order in a, was worked out in 



41 



7l-s>2 



2|p|V, 



POO 

/ dp'dk'6i\p\~p' -k'hip;p\k')[Uil + fp,)il + h)- + U)] , 

^0 

(35) 

and is given in terms of the splitting function for g — > gg. In general this splitting function 
involves the solution of an integral equation which includes the LPM effect. However, in the 



^ In fact we find greatly improved convergence of tlie large- momentum behavior, both in terms of basis set 
size and numerical integration precision, by changing the test functions of [35| to a set which show the 
correct large momentum asymptotic behavior by multiplying 4)2. ..n defined in Eq.(2.32) of the reference 
by 
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deep LPM regime |42j where In (p) can be treated as small, the following leading log result 



for the splitting function can be obtained, 

llM (1 - = -(^V3m ^^^^_^^j3/2 • (36) 

The above splitting function contains the transport parameter q, which characterizes the 
typical transverse momentum squared transferred to the particle per unit length. With the 
above splitting function we show in Appendix [C] that the solution of the off-equilibrium 
distribution function at asymptotically large momentum is 

Up)^^^p'^'. (37) 
4- Summary of weakly coupled pure glue QCD 

Let us now summarize some of the main features of the off-equilibrium dynamics of pure 
glue QCD at weak coupling. In the previous three sections we looked at the behavior of 
the off-equilibrium correction for pure glue QCD in various approximation schemes, deriving 
asymptotic behavior in each case. These asymptotics are listed in Table [B In this section 
we wish to focus on the phenomenologically more interesting region where the equilibrating 
parton has intermediate energies {p ~ lOT). In this case one must resort to numerical 
solutions of the Boltzmann equation which we present in Fig. [3l 

To summarize Fig. [3l we will discuss the curves from top to bottom starting with the 
"Quadratic" curve. In the leading log{T/mo) approximation the linearized Boltzmann equa- 
tion simplifies to a differential equation, Eq. (1271) . The numerical solution to this has been 
worked out in jsi, 3^, 38] and is well described for all momenta by the asymptotic quadratic 



form, X = p^/2T/x. The numerical result will be presented in a forthcoming work 32[, and 
for now we show the quadratic result as the solid blue line (color online). Next we consid- 
ered QCD with the 2 — i- 2 gluon scattering matrix element at leading order. The agreement 
between the asymptotics derived in the previous section and the numerical solution can be 
found in Appendix [Bl At intermediate momentum we show the numerical solution of the 
Boltzmann equation without inelastic processes as the data points under the curve labeled 
"Coll." . The solid curve is the result of a power law fit at intermediate momentum, x oc p^'^. 
In the leading order (LO) treatment when bremsstrahlung is included, we find further equili- 
bration of the gluons and our numerical results are reasonably described by the fit, x p^-^s. 
Finally, the linear ansatz is also shown in Fig. [3] for comparison. 

One can now ask how the observed x p^'^^ of pure glue at leading order will affect the 
viscous corrections to elliptic flow. First of all, as we have already shown, the integrated 
V2 will change marginally. The differential V2, on the other hand, will be largely affected at 
higher pt- This result is shown in Fig. H] along with the quadratic and linear Ansdtze for 
comparison. 

The above considerations have shown that the relaxation of the high energy tail of the 
distribution is largely controlled by energy loss. The low / intermediate momentum region 
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2 4 6 8 10 12 14 16 



p/T 

FIG. 3: The points are from the numerical solution of the Boltzmann equation for pure glue at 
leading order (LO) and without 1 -f-)- 2 processes (i.e. Collisional energy loss only). The lines are 
X oc p^^" for a = 2, 1.6, 1.38, 1 going from top to bottom. 




PT[GeV] 

FIG. 4: V2{pt) for a perturbative gluon gas at leading order. The linear and quadratic Ansdtze are 
shown for comparison. Run parameters can be found in Fig. [TJ 

is constrained by the shear viscosity via Eq. f|T^ . The strength of the off equilibrium correc- 
tion is controlled by two non-perturbative parameters: 1] at low momentum and q at high 
momentum. This is clearly seen by looking at the forms of x we have found for pure glue 
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FIG. 5: The curves at lower momentum are xip/T) for a perturbative gluon gas at leading order 
for three values of the non-perturbative parameter r]/s. The curves at higher momentum show 
the asymptotic forms of x fo^ three values of the non-perturbative parameter q. There must be a 
consistency between 7]/s and q in order that the curves merge at intermediate momentum. 



QCD at leading order, 



xip) 



2.84r; -1.38 
sT P 
0-7 „1.5 



5<P<10 



In ^{p) < 1 



(38) 



In Fig. |5]we show plots of x various choices of the non-perturbative parameters 7]/s and 
q/T^. The main point to take away is the need for a consistency between rj and q such that 
the low and high momentum regions of x can merge smoothly into one another. The three 
values of q/T^ = 10, 16, 60 we have chosen reproduce the experimentally observed Raa [43 
when convoluted with the Higher Twist 



44 



39, 40, 46 



energy 



AMY [291, lail, 1451 and ASW 
loss models respectively. It appears to be difficult to reconcile the discontinuity of x between 
the lowest shear viscosity rj/s = 0.08 and smallest value of q used in modeling heavy ion 
collisions. 



D. Hadron gas 

One might also ask what scattering behavior is expected at lower temperatures, in a 
hadron gas. How do the highest energy hadrons equilibrate, as a function of hadron energy? 
A complete study requires understanding the energy-dependent hadron-hadron cross section, 
which has nontrivial energy dependence and must be determined from experiment. However 
we should be able to say something about the high momentum behavior. 

In hadron-hadron scattering, the inelastic branching fraction rises with increasing s, dom- 
inating the cross-section for kinetic energies well above Aq^^. Since generically no daughter 
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in an inelastic collision carries more than half the energy of the initial high p particle, we 
can take scatterings to be momentum randomizing (the relaxation time approximation is 
sensible), especially for the highest energy hadrons. The relaxation time is then controlled 
by the scattering rate, tr ~ ncr, with n the hadron number density and a an averaged total 
hadronic cross- sect ion. So what is the behavior of the total hadronic cross-section? At low 
momenta it is complicated by resonances but at large momenta there is universally a rising 
total cross-section. Therefore the relaxation time tr{E) should naively involve a small or 
zero power of that is, a ~ 1 is expected, at least for the very high energy tail.^ Certainly 
we do not expect a = 0. However any more detailed discussion must be either model or 
data driven and lies outside the scope of this paper. 

IV. MULTI-COMPONENT PLASMAS 

The plasmas just considered are treated as single-component, in the sense that all degrees 
of freedom are related to each other by symmetries (spins by parity, colors by gauge invari- 
ance). The quark-gluon plasma is a multi-component plasma. Treating as small and 
TTic as large, the three light quark types behave the same, but the gluons behave differently 
from the quarks. Similarly, the hadronic plasma present at lower temperatures contains 
both baryons and mesons, each of several types. The different components generically have 
different departures from equilibrium, that is, Xguark 7^ Xgiuon, which would manifest as 
different viscous corrections to their pt spectra. In particular, we will argue that faster 
equilibration for baryons than for mesons can give a simple explanation for the "constituent 
quark scaling" jsl-lsj observed in V2{pt) for mesons and baryons, without invoking any model 
of coalescence. 

A. Quark-Gluon plasma 

We now consider a two component gas of quarks and gluons and label the distribution 
functions with subscripts q and g respectively: 

= -np{l + np)xg{p)p'p' {diUj) , 
^fq{p) = -np{l - np)xq{p)p'p' {diUj) . (39) 

For use in hydrodynamic simulations we will again fit the off-equilibrium component of the 
quarks' and gluons' distribution function to the following power law, 

X,ip) = C,{a,)p^'^^ . (40) 

The results of the numerical solution of the Boltzmann equation for the two component 
case are shown as points in Fig. |6l The solid curves are the results of the fit done at 
intermediate momentum (5 < p < 15) with the result ~ ~ 0.62. 



^ Froissart behavior a cx In^(s) suggests tr cx In ^(p). 
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FIG. 6: Off-equilibrium correction for the case of a perturbative two-flavor QGP evaluated at 
leading order. The sub-figure shows the ratio of the quark to gluon correction which asymptotically 

approaches Xquark / Xgluon ~ 1-7. 

In order to solve for the two constants {Cq and Cg) we need two constrains. The first 
constraint relates the coefficients Cq^g to the shear viscosity, 



The sum is over quarks and gluons with degeneracies Vg = 2d a = 16 and Vq = ^dfNf = 24. 
The second constraint comes from fixing the ratio of Xq/Xg to the numerical solution of the 
Boltzmann equation. This ratio is shown in Fig. [6] and at large enough momentum {p > 5) 
we find 



The explicit computation of the two coefficients (C^ g) in terms of the above ratio and rj/s 
is worked out in Appendix |Dl 

In Fig. [7] we show the elliptic flow of quarks and gluons. Note the larger suppression 
for quarks as the gluons are forced into equilibrium much quicker. This quicker relaxation 
can not simply be explained by naively assuming Casimir scaling, Xg/Xg ~ Ca/Cf = 2.25. 
Instead this ratio involves a playoff between the faster equilibration rate of gluons and the 
tendency of identity changing processes qq -H- gg, q -H- qg, g ^ qq to equilibrate disequilib- 
rium between the quarks and gluons. This ratio is evaluated analytically at asymptotically 
large momentum in Appendix O 

The distinct quark and gluon elliptic flow is completely due to the different viscous 
corrections, which in turn is related to the different relaxation rates of quarks and gluons. 




(41) 



Xq 

Xg 



1.70. 



(42) 
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Pj [GeV] Pj/n [GeV] 

FIG. 7: Left: Elliptic flow of quarks and gluons. Right: Both V2 and px scaled by n=3,2 for gluons 
and quarks respectively. Run parameters can be found in Fig. [TJ 

Let us note that if we scale both the V2 and pt of gluons by three and quarks by two, the 
result is a "universal curve" as shown in the right plot of Fig. [3 The observed scaling 
is completely accidental, but it led us to consider the possibility of finding similar scaling 
behavior in a meson / baryon system due to differences in the relaxation rates. This is 
discussed in detail in the next section. 



B. Two component meson/baryon gas 

The QCD matter immediately before freezeout is certainly not a weakly coupled quark- 
gluon plasma, but it might be described as a hadron (resonance) gas. Just as for the 
quarks and gluons, there is no reason to think that the mesons and baryons should show 
the same efficiency in equilibrating. But rather than claim a specific model for the partially 
equilibrated state of such a system, we will just do some phenomenology to see how different 
thermalization rates could affect the observed species-dependent elliptic flow behavior. To 
study the hydrodynamics of this system we switch from the Ideal gas equation of state to a 



lattice motivated 47(1 . Further details of the simulation are presented in Appendix |Al 

We consider a meson / baryon gas whereby mesons and baryons have the off-equilibrium 
corrections and fb respectively, 

Sfmip) = -ripil + np)xmip)p'p' {diUj) , 

Sfbip) = -np{l - np)xb{p)p'p^ (diUj) . (43) 
We assume both species have the same power-law correction to spectra. 



Xm{p) = Cm{a)p 

XbiP) = Cb{a)p^-^ , (44) 
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FIG. 8: Elliptic flow of Ks mesons and A baryons from viscous hydrodynamics with radiative or 
quadratic Ansdtze. The run parameters are t]/s = 0.16, Tf^zout = 150 MeV and Lattice EoS. Further 
details are in Appendix [Al The data is from the STAR collaboration 0]. 

but we allow for different coefficients (Cm/Cb) which we will choose in order to give reason- 
able agreement with data. For simplicity, we will consider two different Ansdtze: quadratic 
[a = 0) and radiative {a = 0.5), and take the following ratios which, as we will show, fit the 
data rather well: 

Cm f 1.6 quadratic. 



Cb \ 1-4 radiative. ^^^^ 

Finally, the numerical values of the coefficients can be identified with the shear viscosity 
through 

^ = ^ E ^^^rn/, I -J^^f"'<E^) [1 ± ' (46) 

a=TT,K,... ' 

where the sum extends over all mesons/baryons having M < 1.8/2.0 GeV respectively. This 
choice reproduces the lattice parametrization of the equation of state below T = 160 MeV. 
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We find the following values for the coefficients at our freeze-out temperature of T = 150 
MeV, 



Cm 


= 1.053 1 




a 


= 0.658 J 




Cm 


= 2.661 1 




Cb 


= 1.901 J 


\Ci 



quadratic, (47) 

radiative. (48) 

Before computing particle spectra we would like to make an aside about the way elliptic 
flow is computed. By definition V2{pt) is given by 

, , fd(l)cos(2(p)(dN + 5dN) 

where dN is short for dN/[dpTd(f\ and 5dN is the first viscous correction to this. In the 
above expression the viscous correction to the phase space distribution, 5dN , occurs both in 
the numerator as well as in the normalization from the denominator. Since we have restricted 
the viscous correction to be linear in gradients of field quantities we should therefore require 
that V2 be computed to the same order. We therefore expand the denominator 

/ rf0cos(20) dN + 5dN f d^ 6dN f rf0cos(20) dN 
dN iJd^lNf ' 

so the expression retains terms to first order in 6f only. In the following we will show both 
the expanded and unexpanded expressions for V2, shading the region between the two results 
in order to give an estimate for the uncertainty in the gradient expansion. The upper limit 
of the band corresponds to Eq. P9|) while the lower limit is Eq. f l50|) . In figures where the 
uncertainty band is omitted the plotted curve corresponds to Eq. fH9l) . 

Let us now discuss how the different Ansdtze fare with the experimental data. We have 
chosen rj/s = 0.16 in order to give reasonable agreement with the data in the transverse 
momentum range 1 < pr [GeV] < 2. The f2(pr) spectra for Ks and A are presented 
in Fig. [8] using either the radiative (p^^) or quadratic {p^) Ansatz. For < 2 GeV we 
find good agreement between the viscous hydrodynamic results and data. If one included 
hadronic rescattering the low momentum component of the A V2 would be pushed out 
towards higher px giving better agreement with the data. Above 2-3 GeV large differences 
between the radiative and quadratic Ansdtze are realized. We must warn that at higher pt 
one cannot make a direct comparison with data since a larger fraction of the yield will come 
from fragmenting partons, which have not been included. In addition, the hydrodynamic 
description starts to break down at larger p^- Regardless, one must keep in mind that for 
large enough momentum {i.e. Pr ^ 2 — 3 GeV) the two Ansdtze used here are clearly 
discernible and the choice of Ansatz could in principle lead to differences in the extracted 
viscosity. 

We would now like to investigate whether we observe a meson / baryon scaling, similar 
to the accidental quark / gluon scaling we found from first principles earlier. For clarity, we 
again present the above results with mesons and baryons on the same figure. This is shown 
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FIG. 9: Left: V2 of Ks and A. Right: Constituent quark scaling of V2{pt)- Run parameters can 
be found in Fig. [5J The data is from the STAR collaboration [6| and is plotted in (My — m)/n 
as s ugg ested by the PHENIX collaboration [3] ■ (Recent PHENIX data presented at Quark Matter 



48l . |49|] clearly deviate from constituent quark scaling above {Mt — m)/n ~ 1 GeV.) 



for both radiative and quadratic Ansdtze in Fig. |9l The figures show the corresponding 
results with both f 2 and pt re-scaled by the number of constituent quarks. The scaling of 
the data is the well-known phenomenon of constituent quark scaling. We find that viscous 
hydrodynamics reproduces this "universal curve" as well. This is due to the difference in 
relaxation rates between mesons and baryons, which was treated as a free parameter. The 
possible microscopic origin of this ratio is discussed further in Section |Vl 

We should also point out that this relaxation time scaling is fairly robust to changes in 
the equation of state. While changing the equation of state will clearly affect the r/ = 
behavior, these changes will only have modest modifications to the viscous correction to the 
distribution functions. The qualitative feature that species with smaller relaxation times 
have a stronger elliptic flow is borne out by Fig. [7] (a quark gluon plasma equation of state) 
and Fig. [9] (a lattice equation of state). Further study of the equation of state is left to 
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future work. 



V. SUMMARY AND DISCUSSION 

In this work we have presented a systematic study of the first viscous correction to the 
thermal distribution function. All simulations of viscous hydrodynamics so far have used 
the quadratic Ansatz 

X{p) oc , (51) 

but this is only an educated guess. 

First we studied the form of 6f (or x{p)) ^ momentum dependent relaxation time 
approximation and derived a simple formula 

x{p) = rM^. (52) 

Examining this formula we considered two special cases tr oc p (where the equilibration 
time is proportional to energy) and r^j = const (where the equilibration time is independent 
of energy). These give rise to quadratic {x{p) ^ P^) and linear (x(p) p) dependence on 
momentum, as is summarized in Table [H We expect that, provided QCD is describable in 
terms of quasi-particles, the first viscous correction should lie between these cases. Fig. [2] 
compares these two extreme limits for the functional form of the viscous correction. It is 
important to emphasize that the two simulations have precisely the same shear viscosity. 
Comparing our results for the elliptic flow in these two theories, we see that the integrated 
elliptic flow f 2 is largely determined by the shear viscosity, while differential quantities such 
as V2{pt) at high px depend on the equilibration rates at high momentum. The integrated 
elliptic flow is determined to a large extent by the hydrodynamic variables e^u^^n^^. (An 



explicit formula relatin g vo to e, and vr'^'" is given in [27| which in turn was motivated by 
earlier observations (isTEtI IsoI. 51 |- ) 

The quadratic ansatz is valid only for fairly specialized theories. For instance, examining 
Table [T] we see that scalar theories follow this Ansatz. The reason is that the cross-section 
falls as so higher-energy particles see a more transparent medium and equilibrate more 
slowly. 

For different reasons the quadratic Ansatz is also valid in a soft scattering approximation 
to high temperature QCD (see Row 4 of Table [T]). In this limit, which treats log(T/m£)) 
as an expansion parameter, soft gT collisions lead to the momentum diffusion and drag of 
hard gluons. If the momentum diffusion is independent of particle energy and the drag is 
constant, we get the quadratic Ansatz. If the momentum diffusion increases logarithmically 
with particle energy, we find a logarithmic correction to this Ansatz (see Row 5 of Table [T]). 
A formula which summarizes the asymptotic form of both of these cases is 

2 

= wmw, ' ^^^^ 

where {dE/dt)^ is the rate is energy loss of a particle with momentum p (see Eq. fl26l) and 
Eq. (!32|) for explicit formulas in certain limits). 
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Model 


Physics 


Formula 


Relaxation time, tr(x p 
Relaxation time , tr = const 


Relaxation time grows with particle 
momentum. 

Relaxation time independent of momentum. 


x{p) « 
x{p) « P 


Scalar theory 


Randomizing collisions which happen rarely 


x{p) « P^ 


QCD Soft Scatt. 
QCD Hard Scatt. 
QCD Rad. E-loss 


Soft q ~ gT collisions lead to a random walk of 
hard particles. 

Hard q ~ ^JpT collisions lead to a random walk 
of hard particles. 

Radiative energy controls the approach to equi- 
librium. In the LPM regime q controls the ra- 
diation rate. 


x{p) « P^ 

2 

Xiv) OC 1 T—FP\ 

A.Xl' J ^ log(p/T) 



TABLE I: Summary of the functional dependence of the departure from equilibrium on the theory 
and approximation considered. 



However, the effect of bremsstrahlung completely changes this picture. A naive (Bethe- 
Heitler) treatment of radiative energy loss would lead to a relaxation rate independent of 
momentum, but including the LPM effect, the viscous correction behaves asymptotically as 

3/2 

X(p) = 0.7^. (54) 

This formula is summarized in Row 6 of Table [T] and provides a concrete connection be- 
tween viscous corrections and radiative energy loss which is further explored in Fig. |5] and 
surrounding text. 

From a phenomenological perspective, the LPM effect is not entirely dominant and col- 
lisions are important in the relevant momentum range. A phenomenological fit to numer- 
ical results for the first viscous correction, including both collisions and collinear radiation 
without making the strict LPM approximation, shows that the first viscous correction is 
reasonably well described by the following phenomenological form: 

x{p)^Cp^-^^. (55) 

Fig. m compares this functional form to the linear and quadratic Ansdtze motivated by the 
relaxation time approximation. We see that the general expectation from high temperature 
QCD is that in the relevant momentum range the first viscous correction is slightly closer 
to the linear rather than the quadratic ansatz. 

We next studied two component plasma starting with a two component plasma of quarks 
and gluons. Since the relaxation rates of the quarks and gluons are not the same the two 
components do not have the same distribution function. At high momentum an analysis of 
collinear splittings g ^ gg, g ^ qq, q ^ gq shows that both the quark and gluon distribution 
behave as p^/^. However the ratio of the quark and gluon viscous corrections approaches a 
constant 

^^ 1.70. (56) 

Xg 
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The constant is determined by the ratio of Casimirs Ca/Cf = 9/4 and the dynamics of the 
QCD sphtting functions. It also depends weakly on the number of quark flavors and we 
have quoted the two flavor case. 

Motivated by this example, we have postulated that the baryon and meson components of 
the medium have different equilibration rates. Indeed, there is no reason to expect that these 
species would equilibrate at the same rate. Then we fitted (by eye) the ratio of relaxation 
rates to reproduce the baryon and meson elliptic flows. If the ratio of relaxation times is 

^^1.5, (57) 

Xb 

meaning that baryons relax to equilibrium 1.5 times faster than mesons, then the resulting 
viscous hydrodynamic calculation effortlessly reproduces the universal "constituent quark 
scaling" curve. Physically what is happening is that in ideal hydrodynamics the baryons 
and mesons have approximately the same elliptic flow which is approximately described by a 
linear rise in ttit- The viscous correction then dictates that the baryons will follow this ideal 
trend 1.5 times farther than the mesons. Although it is not obvious from the data shown 
in Fig. [H the data do not show scaling above (mr — Mo)/nq ~ 1 GeV, i.e. the last Lambda 



point is a fluctuation upward (This is seen quite clearly in recent PHENIX data[48|, |49|.) It 
is interesting that the data also deviate from hydrodynamic predictions above this point. 

It is tempting to speculate as to the microscopic origin of the factor of 1.5. The baryons 
and mesons in the 2 — 3 GeV region are produced in the complex transition region where 
the energy density decreases from 1.2GeV/fm^ to 0.5GeV/fm'^. In this range, the temper- 
ature decreases by only AT ~ 20 MeV. However, the hydrodynamic simulations evolve this 
complicated region for a significant period of time, r ~ 4fm -H- 6.5 fm, and the hadronic 
currents are built up over this time period. The interactions are probably quite inelastic 
and are not easily classified as hadronic or partonic in nature. The additive quark model 
was used to describe high energy total cross sections which are similarly inelastic \5^. It 
predicts the ratio of high energy nucleon-nucleon to pion-nucleon (as well as pion-nucleon 
to pion-pion) cross sections to be 3/2 in reasonable agreement with the experimental ratio. 
Perhaps similar physics is responsible for the different baryon and meson elliptic flows. In 
fact, the splitting of the baryonic and mesonic elliptic flows was predicted at least qualita- 



tively by UrQMD which implements the additive quark model [53|]. On the other hand, the 
factor of 1.5 in the relative relaxation times could be simply a combination of dynamical 
and group theoretical factors of accidental significance. 

In summary, a species dependent relaxation time provides a coherent and physically trans- 
parent explanation for the complicated trends observed in the elliptic flow data measured 
at RHIC. 
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Appendix A: Details of hydrodynamic description 

The initial condition of the hydrodynamic evolution is set by a Glauber model and the 
energy density is proportional to the number of binary collisions. More specifically we take 

e(ro = 1 fm, x, y) = Ebc x (Al) 

where Ebc = 22.735 is the energy per binary collision and a^N = 40 mb is the inelastic 
nucleon-nucleon cross section. 

In this work we will use the following evolution equation for tt*-', 

= --{tt'' - W'') - ^Ti'^dky!" + 7^ki^UJ% + -7rfc(,4, (A2) 



which is identical to the stress tensor used in [17j. Other possibilities are also possible [5 



which will not change the results of this work on a qualitative level. In the above expression 
uJij = djUi — diUj is the vorticity and = 3r]/{4:p). There is one technical detail that 
warrants discussion. At large transverse distances the viscous pressure tends to become 
larger than the ideal pressure and the equations become unstable. It is therefore necessary 
to cutoff our auxiliary tensor when it becomes large. More precisely we take 



TT 



where TrTr^ = ^/ nf-^ + 1122 + '^'33 and k ^ 0.1/ (ap). 

In the first part of this paper we consider an ideal gas equation of state, p = l/3e. 
For a two flavor QGP the ideal Stefan Boltzmann gas gives e = 12.71T^, which roughly 
corresponds to the e/T^ relation found on the lattice. (For the highest temperatures in the 
simulation it is above this value and for the lowest temperatures in the simulation it is this 
value). We have decided to use the same e/T"^ ratio for both the gluon gas and quark + 
glue simulations in order to get the fairest possible phenomenological estimate for the size 
of the viscous corrections in a realistic heavy ion event. 

For simulations using the ideal gas EoS the freeze-out contour is taken at constant e^^out = 
0.6 GeV/fm^ corresponding to a temperature of 140 MeV. The default impact parameter is 
7.6 fm and the shear viscosity to entropy ratio is rj/s = 0.08. 

In the second part of this paper where we compute spectra of a meson/baryon gas we 
have used a lattice motivated equation of state [471]. In this case the freeze-out surface is 
set by ef.^out = 0.24 GeV/fm^ corresponding to a temperature of 150 MeV. We have used 
a default impact parameter of 6.8 fm corresponding to a centrality class of 10-40% and a 
shear viscosity to entropy ratio of?7/s = 0.16. 
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Appendix B: Collision Integrals 



In this section we will give the details leading to Eq. fl24|) for a scalar theory and Eq. f l3T|) 
for pure glue. 

1. Scalar theory 

Our starting point is Eq. ( 122|) . Substituting the form specified in Eq. (!23|) into this 
equation yields in a Boltzmann approximation 

T p 2p J {2n)^2k 2 J {2ny ^ ^ 

X 6\iP+K) -P'- K') {P,P, + K,K, - P'^Pl - K'^K) , (Bl) 

where {d^u'^) is the Lorentz invariant extension^ of (d^u^). The integrals over P' and K' 
are Lorentz covariant and can be performed by standard tricks; the {P^Pu + K^Ky) term 
can be factored out, leading to /p, ^, = l/Svr, while the integral over P'^Pl + K'^K^ must 
return a rank-2 tensor depending only on [P+K]^. There are only two such tensors, and 
contraction with g^^'^ and {P+K)^^{P+KY establishes that 

/ ^^^^'^(^")'^(^")'^'((^+^) KK+KK) 

= ^ {A{P+K),{P+K), - {P+Kfg,.) . (B2) 

The integral equation becomes 

{diUj)piPj C (diUj) X'^ f d^k 



T 327rT3 



where we used that {d^u'^) is traceless, {d^u") g^^ = 0. Performing the k angular integration 
in the plasma frame, the pikj terms integrate to zero; so does the kikj term, because (diUj) 
is traceless. Performing the trivial radial integration, we find Eq. (^^. 



2. Pure glue 

Our goal here is to derive Eq. ( 13T|) and Eq. ( I32l) . Our starting point is the collision 
integral Eq. with matrix elements given by Eq. f l2^ 



(diUj) = / ripHkil + nfcO(l + np>)Tpk^p>k' [xip] + x{k) - x{p') - x{k')] ■ (B4) 

Jp'kk' 



^ Specifically, defining the projector onto the local rest frame A^" = g^" + u^u"^, we have 
In the local rest frame implicit here, we have {d'^u^') — and (d^u^) — for fj,,!^ = 1..3. 
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Using the definition, x{p) = x{p) ip^'P " <5*V3) (diUj), one can pull out the common factor, 
(diUj). The remaining integral on the right hand side (called P^) must have the form 
pj = [p^pi — 5'^^ /?)) since this is the only symmetric traceless tensor which can be 
constructed out of p and 5*-' . Straightforward analysis then shows that 



pnp{l + 
T 



where for instance 



npnk{l + nk'){l + npi)Tpk^pik' (B5) 

p'kk' 

X [x{p) + xik)P2icos 9pk) - x{p')P2icos9pp>) - A(cos ^pfc/)] , 



P2{cosekp) = ^ (^^fP - U'^^ [k^y - U^^^ , (B6) 

is the second Legendre polynomial. 

We will evaluate this integral in a leading log(p/T) approximation. Asymptotically, the 
momenta p and p' are large, while k and k' are of order the temperature^. In this limit 
we can make the Boltzmann approximation, (1 + rip/) — )■ 1, and can treat p' as close to p. 
Specifically we take 

cos V = 1 + 77^-1, (B7) 
2pp 

and then write 

dx 

X(p) - X(p) - -^t^- (B8) 

We also note that k and k' are close to T and therefore x{k) and x(^') are small. Then we 
can write Eq. ( ]B6I) as 

P _ ^ /^\ 
dp\dt/^' 

where the average energy loss rate for a particle with momentum p is 



(B9) 



dE\ f , , , , 

— ) =-/ Tpk^p>k'nk{l+nk')uj , (BIO) 

/ p J p'kk' 

i. e. the energy loss is the transition rate weighted with the energy transfer. The energy loss 
to leading log(p/T) has been determined by Bjorken 36|| and Braaten and Thoma 3J] and 
reads 

'fb^'-(^)- 



We will discuss the region of phase space where t = — (P' — P)^ is small. Since the particles are identical, 
there is also an equal contribution where u = —{K' — P)^ is small, i.e. when p and k' are large and p' 
and k are of order T. Our original definition of the transition rate includes a 1/2 symmetry factor for 
the identical particle final state. To ease the discussion in this section, we will simply drop the symmetry 
factor and neglect it-channel contribution. 
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FIG. 10: Off-equilibrium correction for the case of collisional energy loss. The points are from 
the numerical solution of the linearized Boltzmann equation and the curve is the asymptotic form, 
Eq. ()B12p . Specifically the curve is a one parameter fit to the Nc = 3 form, (247r/9) p^/\og{p/C), 
with fit parameter C^^ = 1.3. 



One can verify that when terms suppressed by log(p/T) are dropped we have 

'^'f" = LI T{dE/dt), " tWJwM ' 

Fig. [To] shows a fit based on Eq. ( 1B12I) which does a reasonable job in reproducing our Nc = 3 
numerical results at high momentum. 

For completeness we will rederive Eq. fIBlip . In order to evaluate the phase space integrals 
over Tpk~^p/k/ we use the "t-channel parametrization" of 35| . Following the logic that leads 
from (A. 14) to (A.21) of this work, we write the phase space as 



p'kk' (27r)'^16p2 Jpg 



Tpk^p'k' = , / dk I d(t) I dq I dtu\M\\ (B13) 



PS 



where the momentum transfer is q = p' —p = k — k' and the energy transfer is u = p' —p = 
k — k' . The vector q is taken along the z axis and the vector p lies in the z — x plane. The 
angle is the azimuthal angle of k with respect to the z, x plane. The energy transfer and 
momentum transfer are restricted to the available phase space 

0<^<A:, (B14) 



< ^— ^ <p, (B15) 



which is also exhibited in Fig. [TT] 
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Phase Space (PS) 




-p 



FIG. 11: The available phase space for the collision integrals in Eq. ()B13|) . The band shows the 
dominant region of the integration in a leading \og{p/T) approximation. 



Using the definitions of the kinematic variables, the Mandelstam invariants are 

t = -[p' - pf = q^ (B16) 

s = -{P + K^ = ^\{p + p'){k + k') + g2 _ cos0V(4pp' + t){Akk' + t)] , (B17) 
2g^ L J 

u = -t-s. (B18) 



To evaluate the collision integral, we are to substitute these expressions for the Mandelstam 
invariants into the matrix elements and perform the integrals over the phase space. Close 
inspection of the result of this procedure shows how \og{p/T) comes about. First, the 
logarithm comes from integrating over the phase space region where u ~ —q and T <C 
g ^ p as shown by the band shown in Fig. [11] Since the u integral is over the interval, 
—q<u<—q + 2k, the phase-space integral is approximately 



dq du ^ 2k dq . 

PS J PS J^T 



(B19) 



and we may neglect the stimulation factor, (1 + n^/) ~ 1. Second, only the highest powers 
of u and q contribute to the ultraviolet logarithm. The integrated matrix elements with 
these restrictions is 



2n 



-US 





Then the total total transition rate is 
'dE\ _ 1 



[Sg^Cl] / dk2knk 



dq—^{-q) 

T H 



dt I ^ (27r)4l6p2 
Performing the integral over /c, we arrive at the result quoted in Eq. (IBlip 



(B20) 



(B21) 
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Appendix C: Viscous distribution function and q 



In this appendix we derive the form of the viscous distribution function for asymptotically 
large momenta. In doing this we will relate the high px tail of the distribution function with 
the energy loss parameter q. 

The starting point is the Boltzmann equation containing near coUinear splitting processes. 
We neglect 2 ■<-> 2 processes as these will be sub-leading at large momenta. Therefore 



|!^5./a(x,p) = -Cr^[/], (CI) 

hip 



where 

2|p|V„ 



- / dp'«(|p|-p'-A;')7(p;P,^')[/p(l±/p')(l±A')-/p'/fc'(l±/p)] • 

a JO 

(C2) 

In the above expression /„ is the distribution function of species a. The degeneracy factor, 
i^a is 16 for gluons and 6 for quarks. 
Now linearize the collision integral 

^i->2 ^ ^ r dp'dk'6{p -p'- fc')7(p;p', ± ± ^fcO \Xv - Xp' - Xk\ . (C3) 

Jo 

Doing the integral over k' and expanding out the LHS in the typical way we get 

/3np{l±np) — = / dx^{p;xp, (l-a;)p)np(l±na;p)(l±n(i_^)p) [xp - Xxp - Xii-x)p] ■ 

-t^p ^PVa Jo 

(C4) 

Then note at very high momentum 

{l±n^p){l ±n(i_3;)p) 



1 ± Up 



^ 0(1 - x) . (C5) 



Let us now consider a two component plasma of quarks and gluons. Using the Ansatz 
Xq,g{p) = Cq^gp'^~°' we are left with 



P^Ug 



^2^)3 = Ip'"' I' dXllgiP: xp, (1 - X)P) [Cg - CgX'-^ - Cg{l - X)^""] 

V"" dx^^^ip; xp, (1 - x)p) [Cg - C.x^-'^ - C,(l - xf-^\ , 
Jo 



p'^NiUg 



(C6) 



The sphtting functions at leading log order are 



g. n_ ^ ^ - y^^sCEdpNi ^ V(2Cf - Ca) + Cax'' + Ca(1-x)2 [x^ + (l-x)^] 

^qq[Pi ^Pi ^)P) ~ /n_\4 V PQ' 



(27r)4 [x{l-xp' 



q, , , _ y6Q;,Cj.dj.jVf ^^^CA + {2CF-CA)x'' + CA{l-xf[l + {l- 

lqq\Pi-^Pi\^ X)P) — I VP? 



(27r)^ v.. a;[x(l-x)]^/^ 
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The p^/^ behavior here together with the behavior exphcitly on the RHS of Eq. f lC6p 
must cancel the behavior on the LHS of Eq. ( 1C6|) . This fixes 2 = 2 — a + l/2ora = 1/2, 
so xip) c<: p^^"^- This proves the claim in the main text that the asymptotic behavior should 
be a = 1/2. 



For a gluon gas we find 

^3/2^ htt: r, ^^5/2 



= "y^ - - ""l X [1 - x^/^ - (1 - xm dx , 



(C8) 



and 



0.704778 

= — . (C9) 



For a two-flavor quark-gluon gas we find 

_ 0.759158 

asTVq 



1.257913 , 
X, = -—rj^v"^- (CIO) 



The ratio is 



^ = 1.657, (Cll) 

not too different from the ratio 1.7 we found by fitting. This ratio depends on A^f. For 1 
flavor it is 1.702, for 3 flavors it is 1.618, and in the limit of infinite flavors it approaches 
1.128. This diminishing ratio occurs because, at larger A^f, more and more splitting processes 
are g ^ qq and q qg, which equilibrate the numbers of quarks and gluons towards each 
other. 



Appendix D: Two component system 

In this appendix we derive relationships between the off-equilibrium distribution function, 
X and the shear viscosity of a two component system. 

We consider a gas of bosons and fermions as it will have applications to a gas of quarks 
and gluons or a gas of mesons and baryons. 

^ff{p) = -J^f^pil - np)xf{p)p'p' {diUj) , 

^fbip) = -i^bnpil + np)xbip)p'p' {diUj) . (Dl) 
The off-equilibrium correction x takes the form 

Xbip) = C,{T)p'-'^^ , 

Xfip) = Cf{T)p'~^f . (D2) 
The goal is to find values of the coefficients Cb and C/ as a function of T and rj/s. 
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First we define the partial viscosity of eacli species 



Uf r (Pp 



Vf = — I 77^PXf{p)np [1 - Up] , 



15 J {27i) 



'nb = — I 77TZ\^PXb{p)np [1 + Up] , (D3) 



wliicli will yield a total viscosity of 

V = Vf + Vb- (D4) 

For massive particles the phase space integrals must be done numerically, but for massless 
particles, integrating Eq. (]D3P yields 



Cf{T) 



6s/r3-°/r(6 - a/)C-(5 - a/) 



^^^^^ - 35,r3-".r(6 - a,)C(5 - a,) ' 

with C-i^) = = (1 ~ 2-'^"^X(^)- Let us define TZ as the ratio of the partial 

viscosities, 

^ = 7^. (D6) 



Making use of the relations 



Vb 



V = Vf + Vb = il + '^)Vb, 

77r2 271^ 



s 



we find 



^^^^^ ' 1 + ^ j s 6T3-/r(6 - a^)C-(5 - ' 
^'^^^ = ' TT^) I 3T3-.r(6 X,)C(5 - a,) • ^^^^ 
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